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1 .  INTRODUCTION 

Thiebaux  and  Pedder  (1987)  have  categorized 
types  of  objective  analysis  used  in  the  atmospheric 
sciences  into  three  classes:  surface  fitting, 
empirical  linear  interpolation,  and  statistical 
objective  analysis.  The  recursive  filter  method 
described  here  falls  in  the  second  lot,  and  is  thus 
something  of  a  throwback  from  the  currently 
fashionable  third  category.  However,  the  recursive 
filter  is  faster  and  more  economical  and  thus  more 
practical  for  applications  with  small  computers  and 
limited  time.  For  this  reason,  it  has  been 
Implemented  on  the  University  of  Wisconsin  McIDAS  to 
provide  gridpoint  analyses  from  the  real  time  data 
bases  available  on  that  system.  Why  another  empirical 
fitting  technique  when  other  excellent  methods  exist 
(e.g.,  Cressman,  1959  or  Barnes,  1973)?  We  believe 
that  the  recursive  filter  technique  provides  rather 
unique  capabilities  for  handling  domain  boundaries  and 
background  fields,  and  for  locally  variable  scaling 
dictated  by  observation  quantity  and  quality.  Also, 
it  provides  gridpoint  estimates  of  the  reliability 
of  the  final  analysis  as  a  by-product  (somewhat  like 
the  more  elegant  statistical  schemes).  This  is  useful 
for  both  qualitative  and  quantitative  application  of 
the  analyses. 

The  basic  algorithm  has  been  discussed 
previously  in  Hayden  and  Purser  (1989)  with  a  number 
of  examples  to  elucidate  the  capability  of  the  system. 
In  this  paper,  we  amplify  on  some  special  features  of 
the  two-dimensional  system  and  also  discuss  and 
present  examples  of  the  extension  to  three  dimensions. 

2.  THE  TWO-DIMENSIONAL  RECURSIVE  FILTER  OBJECTIVE 

ANALYSIS  ALGORITHM 

This  section  treats  the  analysis  algorithm  as 
the  sum  of  five  parts.  The  first  part  describes  the 
recursive  filter  which  is  simply  a  technique  for 
smoothing  the  information  at  individual  gridpoints  to 
surrounding  gridpoints.  The  treatment  of  boundary 
points  is  included.  The  second  part  describes  how 
random  data  and  quality  indicators  are  interpolated  to 
regular  gridpoints,  so  that  the  recursive  filter  can 
be  applied.  The  third  part  explains  the  quality 
control  algorithm  and  the  update  of  the  background 
field.  The  fourth  part  discusses  the  use  of  an 
initial  background  field,  and  the  *ifth  part  d^rrlbes 
successive  correction  iterations  which  complete  the 
analysis . 


2.1  The  Recursive  Filter 


The  recursive  filter  is  a  smoother  which  is 
applied  over  a  grid.  The  basic  algorithm  is  described 

by: 

B  -  0  B  ,  +  (l-“)  A  (1) 

n  n-1  n 


According  to  (1),  for  gridpoint  n,  a  new  estimate 
B  is  obtained  from  a  weighted  average  of  the 
previous  value  A  and  the  new  estimate  at  the 
preceding  gridpoint  (n-1).  The  weighting 
parameter  can  be  considered  a  smoothing  parameter 
since  it  controls  the  influence  of  the  preceding 
gridpoints.  This  is  more  easily  seen  when  (1)  is 
written  in  the  recursive  form: 


B  -  (l-°)  (A  +  “A 
n  n  n- 1 


.+  am  A  ) 
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In  practice,  (1)  is  applied  repeatedly  across  the 
rows  and  columns  of  the  grid  in  alternating 
directions  to  provide  symmetry.  The  aggregate 
smoothing  is  equivalent  (Purser  and  McQuigg,  1982) 
to  convolving  the  field  with  a  Gaussian  filter  of 
the  form: 

G(x)  -  e-xs/4T>a  /  ((2r)f2rl))!5  (3) 

where  L  is  the  number  of  passes  of  the  filter  and 
*  is  related  to  the  smoothing  parameter  ° 
according  to 

*  -  1/log  0  ;  **  s  a  /  ()-£>)*  (9) 


The  intrinsic  scale  of  the  filter  is  given  by  the 
standard  deviation  (in  units  of  grid  intervals): 


which  can  be  related  to  a  characteristic  wave¬ 
length  R  if  the  gridlength  is  given  by  *  : 

02  -  (R/5)2  -  2L*2  -  2IP  /(l-°  )2  (6) 

In  using  the  recursive  filter  analysis  it  is 
convenient  to  think  in  terms  of  characteristic 
wavelengths,  and  this  is  the  parameter  specified 
by  the  user.  The  actual  smoothing  parameter0  of 
(1)  is  internally  computed  according  to  lb).  It 
is  interesting  to  note  that  the  accrued  effect  of 
the  multiple  passes  is  equivalent  to  a  single  pass 


& 


of  a  "Barnes"  system  which  employs  a  Gaussian 
filter  directly  (Barnes,  1964;  Koch  et  al., 

1984). 

The  preceding  discussion  implicitly 
assumes  convolutions  extending  over  an  unlimited 
number  of  gridpoints  along  a  line.  In  reality, 
boundaries  intervene.  It  is  possible,  however, 
by  selecting  suitable  starting  conditions  for 
each  application  of  the  elementary  filter  (1),  to 
simulate  the  effect  of  an  infinite  continuation 
of  the  grid,  thereby  mitigating  spurious  boundary 
effects.  Because  the  specification  becomes 
increasingly  complex  with  successive  passes, 
boundary  conditions  are  defined  in  the  algorithm 
for  only  the  zeroth  through  second  passes,  after 
which  the  formulation  for  the  second  pass  is 
repeated.  The  derivation  is  included  as  Appendix 
A. 

2 . 2  Initial  Grid  Specification 

The  discussion  above  has  treated  the 
recursive  filter  acting  on  a  gridpolnt  field.  In 
practice,  observations  do  not  coincide  with 
gridpoints,  and  it  is  necessary  to  transfer  the 
data  before  applying  the  recursive  filter.  This 
is  accomplished  using  the  adjoint  of  simple 
bi-linear  interpolation  to  distribute  each 
observation  to  the  surrounding  gridpoints  (see 
Fig.  1).  Note  that  RX+RX1  -  RY+RY 1  -  1. 

A  +  B 

RY 

~  RX  +  0  -  R1X  * 


A  +  R1X  *  R1Y  *  0 


B1  '  B  +  RX  *  R1Y  *  0 

C1  -  C  +  F,1X  *  RY  *  0 

D1  -  D  +  RX  *  RY  *0 

Fig.  1.  Bilinear  weighting  to  distribute 
observations  (0)  to  gridpoints  (A,B,C,D). 

Two  grids  are  created.  The  first 
represents  the  observation  weight  which  is  a 
function  of  both  number  and  quality  (to  be 
described  below)  of  data  influencing  each 
gridpolnt.  The  second  represents  the  difference 
between  the  observation(s)  and  some  background 
estimate  (bi-linearly  interpolated  to  the 
observation),  multiplied  by  the  observation 
weight.  These  grids  are  cumulative  in  the  sense 
that  the  contributions  of  multiple  observations 
within  a  grid  square  are  added  at  the  surrounding 
gridpoints.  A  multiple  pass  recursive  filter  is 
applied  to  both  grids  with  smoothing  specified  as 
described  previously. 

2.3  Quality  Control  and  Background  Update 


Observation  quality  is  assessed  according 
to  the  function 

Q  -  1 / < (1*(D/T) 4) ;  D  -  |0-A [  (7) 

where  0  is  the  observed  value,  A  the  value 
interpolated  in  the  background  field,  and  T  is  a 
tolerance  which  defaults  to  twice  the  standard 
deviation  of  the  observations,  but  which  can  also 


be  specified  by  the  user.  Note  that  an 
observation  which  exceeds  the  tolerance  is  not 
"rejected",  but  only  down-weighted.  OW“Q*W, 
where  W  is  a  relative  reliability  assigned  to  the 
observation,  is  the  quantity  which  is  distributed 
to  and  accumulated  at  the  gridpoints  to  form  the 
grid  of  observation  weight. 

The  analysis  background  field  is  updated 
after  the  accumulated  QW  and  QW*(0-A)  grids  which 
have  been  operated  on  by  the  recursive  filter. 

The  new  value  at  a  gridpolnt  is  given  by 

A1  -  A  +  C*QW(0-A)/C*QV  (8) 

where  C*  represents  the  application  of  the 
multiple  pass  recursive  filter. 

The  field  of  the  denominator  of  (8)  is  the 
grid  of  analysis  quality  weights  mentioned  in  the 
introduction.  Note  that  this  has  the  implicit 
scaling  of  the  recursive  filter  combined  with  the 
implicit  quality  of  the  observations.  The  weight 
at  any  gridpolnt  represents  a  local  scaling  (in 
two  dimensions)  which  may  be  written  as 

R„  -  1/(C*QW) J*  (9) 

n  n 

2.4  Initial  Background 

The  background  field  to  initiate  the 
analysis  is  a  uniform  grid  where  each  gridpolnt  is 
assigned  the  average  value  of  the  observations. 
Thus,  there  is  no  requirement  for  a  prior  analysis 
or  forecast  of  the  analysis  variable.  There  is, 
however,  provision  for  a  prior  estimate,  if  it  is 
available;  but  it  is  not  used  in  the  traditional 
manner,  as  an  initial  grid.  Rather,  uniformly 
spaced  "data"  are  taken  from  the  analysis/ 
forecast,  assigned  a  reliability,  and  treated 
exactly  as  other  observations.  The  density  of 
these  "data"  can  be  selected  by  the  user,  but 
defaults  to  a  value  for  every  fourth  gridpolnt. 
There  is  no  unique  advantage  to  this  over  the 
traditional  approach,  but  it  is  simple  way  of 
insuring  that  both  observational  and  background 
information  are  treated  analogously.  The 
reliability  of  the  background  data  may  be  assigned 
by  the  user,  but  defaults  to  a  value  such  that  the 
sum  of  the  background  data  reliability  is  equal 
to  the  sum  of  the  observational  data  reliability. 

2.5  Successive  Approximation 

The  previous  paragraphs  have  defined  the 
basic  elements  of  an  analysis  where  the  influence 
of  data  is  spread  according  to  a  characteristic 
wavelength  R,  with  a  quality  control  related  to  a 
tolerance  T.  The  complete  recursive  analysis 
algorithm  is  comprised  of  a  number  of  iterations, 
each  using  the  result  of  the  previous  iteration  as 
the  background.  On  each  iteration,  the  character¬ 
istic  scaling  is  changed.  It  would  be  possible  to 
vary  the  scaling  isotropically  by  reducing  R  from 
some  initial  R0  to  successively  smaller  values. 

In  this  way,  R  would  play  much  the  same  role  as 
the  "radius  of  influence"  of  other  successive 
correction  schemes  (e.g.,  Cressman,  1959). 

However,  (9)  allows  the  option  of  varying  R 
locally  to  reflect  both  observation  density  and 
quality  as  determined  from  the  previous  iteration. 
In  practice,  the  two  concepts  are  combined  such 
that  (9)  is  used,  except  that  it  is  bounded  to  be 
r.c  larger  than  a  predetermined,  successively 
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shrinking  maximum  and  no  smaller  than  a  predeter¬ 
mined  minimum  (which  must  be  at  least  two  grid- 
lengths).  Thus,  observation  density,  quality,  and 
a  nominal  "radius  of  influence"  share  in  the  local 
smoothing  of  the  observations.  In  this  way, 
erroneous  data  are  controlled,  and  there  is  no 
attempt  to  analyze  scales  not  resolved  by  the 
local  observation  density.  Details  of  how  the 
scaling  is  bounded,  and  how  the  quality  control  is 
varied  on  successive  iterations  are  given  in 
Appendix  B. 

3.  THE  THREE  DIMENSIONAL  EXTENSION 

The  extension  to  three  dimensions  is  fairly 
straightforward.  On  a  single  smoothing  pass,  the 
recursive  filter  is  operated  backwards  and 
forwards  along  rows  and  columns  as  discussed 
above,  but  at  a  number  of  vertical  levels.  The 
filter  is  then  passed  up  and  down  in  the  vertical 
dimension  at  each  horizontal  gridpoint.  This 
procedure  is  repeated  for  the  specified  number  of 
smoothing  passes.  The  maximum  horizontal 
characteristic  wavelength  R0  is  varied  for  each 
level  since  the  atmosphere  tends  to  show  smaller 
scales  at  lower  levels.  The  horizontal  smoothing 
at  each  level  is  related  to  the  vertical 
smoothing  through  a  vector  (ZMET)  which  may  be 
defined  by  the  operator,  but  which  defaults  to  a 
value  such  that  the  vertical  scale  is  about  a 
hundredth  of  the  horizontal  scale.  Interpolation 
to  the  observation  location  along  the  vertical 
currently  uses  a  six  point  (where  possible) 
Lagrangian  interpolation,  since  this  appears  to 
improve  the  convergence  of  the  analysis  as 
compared  to  the  one  dimensional  bi-llnear 
relationship  shown  in  Fig.  1.  However,  the 
observation  corrections  are  distributed  to  the 
surrounding  levels  using  the  bi-linear 
formulation,  since  experiments  using  the  adjoint 
of  the  Lagrangian  interpolation  have  not  been 
satisfactory,  perhaps  because  of  the  varying 
signs  of  the  coefficients.  He  do  not  feel  that 
this  matter  is  completely  resolved.  If  the 
horizontal  versus  vertical  scaling  were  correctly 
specified,  it  seems  that  the  simple  bi-linear 
relationships  should  work  in  all  dimensions. 

A  background  field  may  be  applied  in  the 
three  dimensional  analysis  by  the  same  technique 
described  above.  Background  "data"  are  taken 
from  whatever  analyses/forecasts  are  available. 

In  the  three  dimensional  case,  there  is  a  clear 
advantage  to  this  technique  over  the  traditional 
method  since  analyses/forecasts  do  not  have  to  be 
available  for  all  vertical  levels,  nor  do  they 
need  to  be  in  the  same  vertical  dimension  to  be 
readily  used. 

4.  EXAMPLES 

Figures  2  and  3  present  examples  of  the 
three-dimensional  analysis  results  at  500  mb  for 
the  situation  of  00  UT,  28  January  1986.  Figure 
2  shows  the  geopotential  derived  for  analysis 
iteration  1,  4,  and  8.  Figure  3  shows  the 
corresponding  (normalized)  reliability  field. 

The  analysis  was  done  at  ten  pressure  levels  with 
a  horizontal  gridlength  of  111  km.  A  prior 
analysis  was  used  to  provide  background  data  of 
uniform  reliability  at  every  fourth  gridpoint. 

The  figures  show  the  "successive  approximation" 
feature  of  the  system,  especially  with  respect  to 
how  the  reliability  field  encloses  the  observa¬ 


tions.  Note  that  the  reliability  field  becomes 
uniform  over  the  Atlantic  as  the  iterations 
increase,  reflecting  that  the  analysis  is  a 
product  of  only  background  data.  A  slight 
fall-off  occurs  at  the  edge  of  the  grid,  and  more 
at  the  corner,  reflecting  that  the  boundaries 
prohibit  uniform,  isotropic  data  coverage. 

Figures  2  and  3  can  be  compared  with  Fig.  4  which 
represents  the  eight  iteration  analysis  without 
background  data.  The  reliability  field  shows 
more  concentration  around  the  radiosonde  stations 
and  a  much  sharper  fall-off.  The  analysis  in  the 
southeast  over  the  ocean  is  designated  as 
completely  unreliable,  and  the  geopotential  here 
versus  that  shown  in  Fig.  2  suggests  that  this 
appraisal  is  accurate. 


Fig.  2.  500  mb  geopotential,  00VT  28  January  1986 

top  to  bottom:  one,  four  and  eight  analysis 
iterations . 
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Fig.  4.  Geopotential  and  reliability  analyses  of 
Fig.  2  without  background  estimate. 


Fig.  3.  Reliability  fields  associated  with 
analyses  of  Fig.  2. 

Figures  5  and  6  represent  a  temperature 
analysis  at  750  mb.  Since  this  is  not  a  standard 
reporting  level,  the  analysis  is  a  product  of 
vertical  interpolation  and  significant  level 
reports.  Temperature  contours  are  shown  in  Fig. 

5  for  an  analysis  which  used  the  significant 
level  data.  The  field  of  differences  from  an 
analysis  which  did  not  use  the  significant  level 
data  is  also  shown.  It  can  be  seen  that  the 
significant  level  data  has  sharpened  the  analy¬ 
sis,  making  changes  of  greater  than  2  K  in  the 
vicinity  of  Florida.  Figure  6  shows  the 
reliability  fields  associated  with  the  two 
analyses.  The  one  without  significant  level  data 
shows  that  the  result  is  almost  entirely  composed 
of  the  background  data,  though  there  is  a  slight 
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Fig.  5.  Top:  750  mb  temperature  analysis  using 
significant  level  data.  Bottom:  difference  from 
analysis  without  significant  level  data. 
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Fig.  6.  Reliability  fields  of  750  mb  analysis 
without  (top)  and  with  (bottom)  significant  level 
data . 

skewing  over  the  land.  One  can  also  see  that  the 
reliability  field  with  significant  data  is  lower 
over  the  data  sparse  Atlantic  than  for  the  500  mb 
geopotential  analysis  (Fig.  3).  This  is  because 
the  background  data  reliability  is  adjusted 
according  to  the  density  of  observations  as 
explained  above. 

5.  SUMMARY 

This  paper  has  described  the  features  of  a 
three-dimensional  objective  method  for  deriving 
meteorological  analyses.  While  the  analysis  is 
not  optimal  in  any  rigorous  statistical  sense,  it 
is  highly  versatile  and  efficient  computational¬ 
ly.  For  example,  the  analysis  described  in 
section  4  comprised  961  horizontal  gridpoints  at 
ten  levels  with  1468  observations.  Using  eight 
analysis  iterations,  it  was  accomplished  in  15 
seconds  on  the  McIDAS  (an  IBM  4381).  The 
three-dimensional  form  of  the  analysis  allows 
observations  at  significant  levels  to  be  easily 
incorporated  (especially  Important  for  cloud 
drift  winds).  The  user  is  able  to  specify  the 
relative  weight  of  observations  and  the 
background  field  independently  at  each  level.  A 
simple  quality  control  procedure  allows  the  gross 
error  of  observations  to  be  identified  and 
rendered  ineffective  by  severely  down-weighting 
these  observations  before  they  can  impact  the 
analysis.  Apart  from  efficiency,  the  main 
strength  of  the  scheme  is  its  ability  to  handle 
the  presence  of  data-dense  and  sparse  areas  in 
the  same  domain,  producing  fields  of  reliability 
which  are  useful  in  the  application  of  the 
analyses . 


6.  APPENDIX  A 

The  boundary  problem  arises  in  specifying 
the  updated  value  of  the  leftmost  gridpoint  for  a 
right-moving  filter,  and  vice  versa,  for  multiple 
passes,  since  there  is  an  implied  "tail"  extending 
beyond  the  grid  from  preceding  passes.  For  the 
zeroth  pass,  no  special  treatment  is  required 
since  there  is  no  "tail"  and  the  initial  adjust¬ 
ment  is  just 

B j  -  (1-c)  A j  (Al) 

On  succeeding  passes,  the  tail  can  be  treated 
explicitly.  However,  the  situation  becomes 
increasingly  complicated,  and  as  mentioned  in  the 
main  text,  the  analysis  algorithm  specifies 
starting  conditions  for  only  the  first  and  second 
passes,  after  which  the  relationship  for  the 
second  pass  is  repeated. 


In  general,  the  starting  conditions  for  a 
right-moving  filter  can  be  derived  by  considering 
the  aggregate  of  multiple  applications  of  a 
left-moving  filter  which  is  expressed  by: 


4  Cl 


“2  An+2 


P  n+p 


P 

+  (1  -  Z  a  .)  A°  (A2) 

i-1  1  n 


where  n  is  a  gridpoint  number,  p  is  the  pass 
number,  the0i  are  some  combination  of  the  basic 
smoothing  parameter  a ,  and  the  A  are  the  grid- 
point  values.  Note  that  for  griSpolnt  0  (to  the 
left  of  the  row)  (A2)  reduces  to 

p  p  p  p 

A  “oi  Ai  *1,  Ai  ♦  ■"  A  (A3) 

O  1  1  Z  L  p  P 

since  the  initial  value  A°  is  not  defined  outside 
the  grid.  A  is  precisely  what  is  required  to 
begin  the  right-moving  filter  at  gridpoint  1. 

Thus,  the  problem  reduces  to  specifying  thea^. 

For  p*l,  (A2)  reduces  to 

‘“I  Cl  +  (1-1>  An  (A4) 


and  thus  o  .  "ex.  For  p*2,  there  are  the  basic 
relationships : 
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Eliminating  A^  and  A^+J: 
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which  is  of  the  form  (A2)  where  a  ,  *  ^  and  a ,  * 
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So  to  begin  the  right-moving  filter  for  p»l 
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Eliminating  Aq  and  A2 
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For  the  left-moving  filter,  the  starting  condi¬ 
tions  are  symmetric  since  the  filters  commute. 


7. 


APPENDIX  B 


The  concept  of  locally  varying  scaling 
which  is  determined  by  observation  quantity  and 
quality  has  been  developed  in  the  main  text.  For 
each  gridpoint,  after  the  first  Iteration  of  the 
analysis,  there  is  a  weight  W  which  is  implicit¬ 
ly  related  to  a  scale  R  ,  depending  on  the 
scaling  of  the  recursive  filter  applied  in  the 
previous  iteration 


1 


(Bl) 


(Note:  R  3  is  used  in  the  three-dimensional 
version.)  R  is  used  to  determine  the  local 
scaling  of  tRe  recursive  filter,  subject  to 
bounding  limits  which  will  be  explained  here. 
First,  however,  it  is  necessary  to  express  the 
way  in  which  the  scaling  is  formulated  in  the 
recursive  filter  algorithm.  As  presented  in  the 
main  text,  the  characteristic  wavelength  of  the 
recursive  filter  is  given  by 


(R/*5  ) 2 


2D*  /(  1-a  )  2 


(B2) 


where  L  is  the  number  of  passes  and  $  is  the 
gridlength.  If  the  basic  recursive  algorithm  is 
reformulated  such  that  ®  «  l-°  , 

B  -  (1-fi )  B  ,  +  B  a 
n  n-1  n 

then  (B2)  leads  to 


where 


6  -  -e  +  ^  e(e+2) 


-e  ♦  (e(e+2)) 


e  -  if  a/Ra  -  if  aW 


(B3) 


(B4) 


For  the  first  iu ration  of  the  analysis,  the  R  of 
(B4)  is  R0  which  is  chosen  by  the  operator  as  the 
largest  characteristic  scale  for  the  analysis. 

For  the  second  and  subsequent  iterations  R  is 
replaced  by  R^*R  where  R  is  given  by  (Bl)  and  R^ 
is  a  tuning  consiant  whlcR  will  be  explained 
below.  The  value  is  constrained  such  that 


R  £ 
m 


R  *R 
f  n 


R0 


(B5) 


RL  +  R  (R0  -  RL) 


where  RL  (chosen  by  the  operator)  represents  a 
limiting  wavelength,  and  R  controls  how  rapidly 
the  characteristic  wavelen|th  is  diminished  on 
each  succeeding  iteration  i  (a  factor  of  .33  is 
the  default  value).  Equation  B5  is  shown  schema¬ 
tically  in  Fig.  Bl.  Note  from  (B5)  that  the  scale 
for  the  entire  analysis  is  essentially  limited 
between  R0  and  Rl. 
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Fig.  Bl. 
limit. 


2  3  4  B  • 

ITERATION  INDEX 

Schematic  of  the  assignment  for  R 


The  factor  R{  is  included  to  give  the 
operator  some  control  over  the  scaling  beyond  that 
offered  by  choosing  the  upper  and  lower  limits  of 
R.  R,  effectively  sharpens  or  spreads  the  scaling 
according  to  whether  it  is  less  than  or  greater 
than  one.  For  noisy  data,  a  value  of  2  is  perhaps 
appropriate,  whereas  for  very  clean  data,  a  value 
of  0.5  might  be  used. 
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